* This do-file creates regression tables and figures for the input experiment
* Paper Tables: 3, 5, S6, S7, S8, S9-S34 (not well-numbered)
* Paper Figures: 6, S3, S9-S13
* NOTE: This will create EITHER Figure 6 (if pic==cm) OR Figure S3 (if pic==sd)
* NOTE: To create the other one, run this file again with main=1, pics=1 and other locals = 0
* NOTE: Also displays PCA construction details, but Table S5 is hand-constructed

* Select which pieces to run
	* Tables: 3, 5, S6, S7, S8
	local main = 1
	* Figure 6 OR Figure S3 (set `pic' below)
	local pics = 1
	* Any of the options below
	local appx = 1
		* Double post lasso table (NOT in paper)
		local dlasso = 0
		* Table versions of Figure S9-S13 (NOT in paper)
		local robust = 0
		* Figures S9-S13
		local rob_pic = 1
		* Tables S9-S34
		local hetero = 1
	
* Choose how to rescale for Figure 6 / S3 (cm or sd)
local pic cm
	if 		("`pic'"=="cm") {
		local scl "mean"
		local fno "6"
	}
	else if ("`pic'"=="sd") {
		local scl "sd"
		local fno "S3"
	}

* Name a bunch of stuff
global B_app B_y1_e B_y2_e B_y3_e B_y1_s B_y2_s B_y3_s B_y2 B_y3
global B_het B_y1_t B_y1_t_i B_y2_t B_y2_t_i B_y3_t B_y3_t_i int_i B_y2 B_y2_i B_y3 B_y3_i

// This stuff should be local, but this makes it easier to run sections of code
global grset graphregion(color(white)) ylabel(, angle(0)) grid(glcolor(gs0) glpattern(dash) glwidth(thin) nogextend) plotregion(margin(0) lcolor(black)) 
global xset transform(* = min(max(@,-1),2)) xscale(range(-1(1)2)) xlabel(0 "0" 1 "Control Mean") 

global estset replace se tex label  mlabels(none) nonotes nonumbers starlevels(* 0.1 ** 0.05 *** 0.01)
global mainstats stats(ctrl_mean r2 N, fmt(2 2 0) labels("Control Mean" "R-Squared" "Observations"))  
global hetstats stats(ctrl_mean het_mean r2 N, labels("Ctrl. Mean ($ I=0$)" "Ctrl. Mean ($ I=1$)" "R-Squared" "Observations"))
global mgset prefix(\multicolumn{@span}{c}{) suffix(}) span
global ecell collab(none) gaps cells(b(fmt(3) pvalue(q_val) star) se(fmt(2) par) q_val(fmt(3) par({[ ]})))

tempfile treats
graph drop _all
snapshot erase _all

*****
* Household treatment assignemnts
*****

use "$out_data/survey_roster.dta", clear
merge m:1 village using "$admin_data/treatment_assignment.dta", keepusing(treat_input treat_ext) assert(match using)
drop if _m==2
drop _m
merge m:1 hhid using "$out_data/survey_drops.dta", assert(match) nogen

* Split land and wealth into above/below median for heterogeneity analysis
foreach v in land wealth {
	quietly sum bl_`v'_win5, d
	gen byte bl_`v'_high = (bl_`v'_win5 > `r(p50)') if !missing(bl_`v'_win5)
}

// Make year*treat interactions
expand 3
bys hhid: gen byte year = _n
* Dummy out treatment interactions
gen byte B_y2 = year==2
label var B_y2 "Year 2"
gen byte B_y3 = year==3
label var B_y3 "Year 3"
gen byte B_y1_t = (year==1)*(treat_input==1)
label var B_y1_t "Treat Yr. 1"
gen byte B_y2_t = (year==2)*(treat_input==1)
label var B_y2_t "Treat Yr. 2"
gen byte B_y3_t = (year==3)*(treat_input==1)
label var B_y3_t "Treat Yr. 3"

gen byte B_y1_s = (year==1)*(treat_ext==1)
label var B_y1_s "Subs. Yr. 1"
gen byte B_y2_s = (year==2)*(treat_ext==1)
label var B_y2_s "Subs. Yr. 2"
gen byte B_y3_s = (year==3)*(treat_ext==1)
label var B_y3_s "Subs. Yr. 3"
gen byte B_y1_e = (year==1)*(treat_ext==2)
label var B_y1_e "Ext. Yr. 1"
gen byte B_y2_e = (year==2)*(treat_ext==2)
label var B_y2_e "Ext. Yr. 2"
gen byte B_y3_e = (year==3)*(treat_ext==2)
label var B_y3_e "Ext. Yr. 3"

save `treats'

*****
* Adoption and Area; Production; Yield
*****
* Kharif adoption (dummy), Kharif area, Rabi adoption (dummy), Rabi area, Zaid adoption, Zaid area
use "$out_data/outcome_adoption_production.dta", clear
collapse (sum) pulse_prod_win5 pulse_area, by(hhid year season) fast

* Make adoption dummy
gen byte pulse_adopt = pulse_area>0 if !missing(pulse_area)

* Make pulse yield (total)
gen pulse_yield = pulse_prod_win5 / pulse_area

* Reshape over seasons and merge in pulse stocks
reshape wide pulse_adopt pulse_area pulse_prod_win5 pulse_yield, i(hhid year) j(season)
merge 1:1 hhid year using "$out_data/outcome_stocks.dta", assert(match) nogen keepusing(pulse_time*)
egen pulse_time_tot = rowmax(pulse_time_*)

* Offset pulse stocks by a year
sort hhid year
by hhid (year): gen pulse_time_prev = pulse_time_tot[_n+1]

* Add in treatment data
merge 1:1 hhid year using `treats', assert(match) nogen


*****
* Regressions: Adoption
*****

*** Regressions: Adoption and Area
local v_adopt pulse_adopt1 pulse_area1 pulse_adopt2 pulse_area2 pulse_adopt3 pulse_area3
	label var pulse_adopt1 "Kharif Adoption"
	label var pulse_area1 "Kharif Area"
	label var pulse_adopt2 "Rabi Adoption"
	label var pulse_area2 "Rabi Area"
	label var pulse_adopt3 "Zaid Adoption"
	label var pulse_area3 "Zaid Area"
snapshot save, label("s_adopt")


if `main'==1 {
eststo clear
cap mat drop P Q

foreach vbl of varlist `v_adopt' {
	local RV $Rvars
	if (`:list posof "`vbl'" in RV') local samp ever
	else local samp srvy

	* Regressions for table
	eststo `vbl': reg `vbl' $B_main $controls if drop_`samp'==0, vce(cluster village)
	mat tmp = r(table)
	quietly sum `vbl' if year==1 & treat_input==0
	estadd scalar ctrl_mean = `r(mean)'
	
	* Save p-values for qqvalue
	mat P = nullmat(P) \ tmp[4,1..3]'
}

* Compute q-values
preserve 
clear
svmat P
qqvalue P1, qvalue(qval) method(yekutieli)
mkmat qval, mat(Q)
restore

* Save q-values to estimation output
local i = 1
foreach vbl of varlist `v_adopt' {
	if strpos("`vbl'","pca") continue
	
	local j = `i'+2
	mat q_val = Q[`i'..`j',1]'
	mat colnames q_val = $B_coef
	estadd matrix q_val : `vbl'
	local i = `i'+3
}

* Output results to table
local col1 "  & Adoption & Area & Adoption & Area & Adoption & Area \\"
local col2 " & (1) & (2) & (3) & (4) & (5) & (6) \\ \midrule "
local rules \cmidrule(lr){@span}
local est_adopt mgroups("Kharif" "Rabi" "Zaid", pattern(1 0 1 0 1 0) $mgset erepeat( `rules') ) posthead(`col1' `col2')
esttab using "$out_tables/Table_3.tex", keep($B_main) $estset $mainstats $ecell `est_adopt'
}

*****
* Regressions: Yield
*****

eststo clear
estimates clear
cap mat drop P Q

local vbls pulse_yield1 pulse_yield2 pulse_yield3
/// NO PCA: only 43 obs have yield in each season
// pca `vbls'
// predict pca_yield, score
// local v_yield "`vbls' pca_yield"
local v_yield "`vbls'"
	label var pulse_yield1 "Kharif Yield"
	label var pulse_yield2 "Rabi Yield"
	label var pulse_yield3 "Zaid Yield"
snapshot save, label("s_yield")

if `main'==1 {
foreach vbl of varlist `v_yield' {
	local RV $Rvars
	if (`:list posof "`vbl'" in RV') local samp ever
	else local samp srvy
	
	* Regressions for table
	eststo `vbl': reg `vbl' $B_main $controls if drop_`samp'==0, vce(cluster village)
	mat tmp = r(table)
	sum `vbl' if year==1 & treat_input==0
	estadd scalar ctrl_mean = `r(mean)'
	
	* Save p-values for qqvalue
	mat P = nullmat(P) \ tmp[4,1..3]'
	
}

* Compute q-values
preserve 
clear
svmat P
qqvalue P1, qvalue(qval) method(yekutieli)
mkmat qval, mat(Q)
restore

* Save q-values to estimation output
local i = 1
foreach vbl of varlist `v_yield' {
	if strpos("`vbl'","pca") continue
	
	local j = `i'+2
	mat q_val = Q[`i'..`j',1]'
	mat colnames q_val = $B_coef
	estadd matrix q_val : `vbl'
	local i = `i'+3
}

* Output results to table
local col1 "  & Kharif & Rabi & Zaid \\"
local col2 " & (1) & (2) & (3) \\ \midrule "
local rules \cmidrule(lr){2-4}
local est_yield mgroups("Yield (Kg./Acre)", pattern(1 0 0) $mgset end( `rules')) posthead(`col1' `col2')
esttab using "$out_tables/Table_5.tex", keep($B_main) $estset $mainstats $ecell `est_yield'
}
	
*****
* Regressions: Production and Stocks
*****

eststo clear
estimates clear
cap mat drop P Q

local vbls pulse_prod_win51 pulse_prod_win52 pulse_prod_win53 
pca `vbls', comp(1)
predict pca_prod, score
local v_prod "`vbls' pca_prod pulse_time_prev"
	label var pulse_prod_win51 "Kharif Production"
	label var pulse_prod_win52 "Rabi Production"
	label var pulse_prod_win53 "Zaid Production"
	label var pca_prod "Production Index"
	label var pulse_time_prev "Harvest Stock"
snapshot save, label("s_prod")

if `main'==1 {
foreach vbl of varlist `v_prod' {
	local RV $Rvars
	if (`:list posof "`vbl'" in RV') local samp ever
	else local samp srvy
	
	// This throws off q-value computation; easier to just delete empty cells manually
	if ("`vbl'"=="pulse_time_prev") local regs "B_y1_t B_y2_t B_y2"
	else local regs "$B_main"

	eststo `vbl't: reg `vbl' $B_main $controls if drop_`samp'==0, vce(cluster village)
	mat tmp = r(table)
	quietly sum `vbl' if year==1 & treat_input==0
	estadd scalar ctrl_mean = `r(mean)'
	
	if strpos("`vbl'","pca") continue
	
	* Rescale relative to control mean for figure
	replace `vbl' = `vbl' / `r(`scl')'
	quietly reg `vbl'  $B_main $controls if drop_`samp'==0, vce(cluster village)
	estimates store `vbl'
	
	* Save p-values for qqvalue
	mat P = nullmat(P) \ tmp[4,1..3]'
}

* Compute q-values
preserve 
clear
svmat P
qqvalue P1, qvalue(qval) method(yekutieli)
mkmat qval, mat(Q)
restore

* Save q-values to estimation output
local i = 1
foreach vbl of varlist `v_prod' {
	if strpos("`vbl'","pca") continue
	
	local j = `i'+2
	mat q_val = Q[`i'..`j',1]'
	mat colnames q_val = $B_coef
	estadd matrix q_val = q_val : `vbl't
	local i = `i'+3
}

* Output results to table
local col1 "  & Kharif & Rabi & Zaid & Index & (Months) \\"
local col2 " & (1) & (2) & (3) & (4) & (5) \\ \midrule "
local rules \cmidrule(lr){2-4}
local est_prod mgroups("Harvest Production (Kgs.)" "Prod." "Harvest Stock", pattern(1 0 0 1 1) $mgset end( `rules')) posthead(`col1' `col2')
esttab using "$out_tables/Table_S6.tex", keep($B_main) $estset $mainstats $ecell `est_prod'

* Make graph for picture
if `pics'==1 {
coefplot `vbls' pulse_time_prev ///
	, keep($B_coef ) xline(0) xline(1) $grset $xset name(gprod) ///
	p1(label("Kharif Production") msymbol(o) color(edkblue) ciopts(lcolor(edkblue) lwidth(thin))) ///
	p2(label("Rabi Production") msymbol(d) color(maroon) ciopts(lcolor(maroon) lwidth(thin))) ///
	p3(label("Zaid Production") msymbol(t) color(eltgreen) ciopts(lcolor(eltgreen) lwidth(thin))) ///
	p4(label("Months with Pulses") msymbol(s) color(brown) ciopts(lcolor(brown) lwidth(thin)))
graph export "$out_figures/Parts/fig_coef_prod_`pic'.png", replace
}
}



*****
* Profits
*****
* Regressions with combined treatment
* Profits, Production Revenue, Sales Revenue, Cost, Land Owned
use `treats', clear
merge 1:m hhid year using "$out_data/outcome_profits.dta", assert(match) nogen

*** Main regression
eststo clear
estimates clear
cap mat drop P Q

local vbls earn_revenue_win5 earn_sales_win5 earn_cost_win5 farm_area_acre_win5
pca `vbls', comp(1)
predict pca_earn, score
local v_earn "earn_profit_win5 `vbls' pca_earn"
	label var earn_profit_win5 "Farm Profit"
	label var earn_revenue_win5 "Farm Revenue"
	label var earn_sales_win5 "Crop Sales"
	label var earn_cost_win5 "Farm Cost"
	label var farm_area_acre_win5 "Farm Area"
	label var pca_earn "Profit Index"
snapshot save, label("s_earn")

if `main'==1 {

foreach vbl of varlist `v_earn' {
	local RV $Rvars
	if (`:list posof "`vbl'" in RV') local samp ever
	else local samp srvy
	
	eststo `vbl't: reg `vbl' $B_main $controls if drop_`samp'==0, vce(cluster village)
	mat tmp = r(table)
	quietly sum `vbl' if year==1 & treat_input==0
	estadd scalar ctrl_mean = `r(mean)'
	
	if strpos("`vbl'","pca") continue
	
	* Rescale relative to control mean
	quietly replace `vbl' = `vbl' / `r(`scl')'
	quietly reg `vbl' $B_main $controls if drop_`samp'==0, vce(cluster village)
	estimates store `vbl'
	
	* Save p-values for qqvalue
	mat P = nullmat(P) \ tmp[4,1..3]'
}

* Compute q-values
preserve 
clear
svmat P
qqvalue P1, qvalue(qval) method(yekutieli)
mkmat qval, mat(Q)
restore

* Save q-values to estimation output
local i = 1
foreach vbl of varlist `v_earn' {
	if strpos("`vbl'","pca") continue
	
	local j = `i'+2
	mat q_val = Q[`i'..`j',1]'
	mat colnames q_val = $B_coef
	estadd matrix q_val : `vbl't
	local i = `i'+3
}

* Output results to table
local col1 " & (1) & (2) & (3) & (4) & (5) & (6) \\ \midrule "
local est_earn mgroups("Profit" "Revenue" "Sales" "Cost" "Farm Area" "Profit Index", pattern(1 1 1 1 1 1) $mgset ) posthead(`col1')
esttab using "$out_tables/Table_S7.tex", keep($B_main) $estset $mainstats $ecell `est_earn'

* Make graph for picture
if `pics'==1 {
coefplot earn_profit_win5 `vbls' ///
	, keep($B_coef ) xline(0) xline(1) $grset $xset name(gprofit) legend(cols(3)) ///
	p1(label("Net Profit") msymbol(o) color(edkblue) ciopts(lcolor(edkblue) lwidth(thin))) ///
	p2(label("Production Value") msymbol(d) color(maroon) ciopts(lcolor(maroon) lwidth(thin))) ///
	p3(label("Sales Revenue") msymbol(t) color(eltgreen) ciopts(lcolor(eltgreen) lwidth(thin))) ///
	p4(label("Total Cost") msymbol(s) color(brown) ciopts(lcolor(brown) lwidth(thin))) ///
	p5(label("Area Farmed") msymbol(X) color(gs0) ciopts(lcolor(gs0) lwidth(thin)))
graph export "$out_figures/Parts/fig_coef_profit_`pic'.png", replace
}
}



*****
* Consumption
***** 
* Regressions with combined treatment
* Pulse consumption, Protein consumption, Pulse stocks
use `treats', clear
merge 1:m hhid year using "$out_data/outcome_consumption.dta", assert(match master) nogen
merge m:1 hhid year using "$out_data/outcome_stocks.dta", assert(match) nogen keepusing(stock*)

foreach vbl in kgs kw5 {
	egen stock_`vbl'_tot = rowtotal(stock_`vbl'_*)
}
egen stock_time_tot = rowmax(stock_time_*)

* Main regression
eststo clear
estimates clear
cap mat drop P Q

local vbls stock_kw5_tot stock_time_tot eat_week_pulse_win5_pc eat_day_prot_win5_pc
pca `vbls', comp(1)
pca `vbls', comp(2)
predict pca_eat pca_alt_eat, score
local v_consume "`vbls' eat_day_pfem_win5 pca_eat"
	label var stock_kw5_tot "Pulse Stock (Kgs.)"
	label var stock_time_tot "Pulse Stock (Mos.)"
	label var eat_week_pulse_win5_pc "Consumption"
	label var eat_day_prot_win5_pc "Daily HH Protein"
	label var eat_day_pfem_win5 "Female Protein"
	label var pca_eat "Consumption Index"
snapshot save, label("s_consume")

if `main'==1 {
foreach vbl of varlist `v_consume' {
	local RV $Rvars
	if (`:list posof "`vbl'" in RV') local samp ever
	else local samp srvy
	
	eststo `vbl't: reg `vbl' $B_main $controls if drop_`samp'==0, vce(cluster village)
	mat tmp = r(table)
	if ("`vbl'"=="eat_day_pfem_win5") quietly sum eat_day_pfem_win5 if year==2 & treat_input==0
	else quietly sum `vbl' if year==1 & treat_input==0
	estadd scalar ctrl_mean = `r(mean)'
	
	if strpos("`vbl'","pca") continue
	
	* Rescale relative to control mean
	quietly replace `vbl' = `vbl' / `r(`scl')'
	quietly reg `vbl' $B_main $controls if drop_`samp'==0, vce(cluster village)
	estimates store `vbl'
	
	* Save p-values for qqvalue
	mat P = nullmat(P) \ tmp[4,1..3]'
}

* Compute q-values
preserve 
clear
svmat P
qqvalue P1, qvalue(qval) method(yekutieli)
mkmat qval, mat(Q)
restore

* Save q-values to estimation output
local i = 1
foreach vbl of varlist `v_consume' {
	if strpos("`vbl'","pca") continue
	
	local j = `i'+2
	mat q_val = Q[`i'..`j',1]'
	mat colnames q_val = $B_coef
	estadd matrix q_val : `vbl't
	local i = `i'+3
}

* Output results to table	
local col1 "  & Kgs. & Months & (Kg./Week) & HH/capita & Female & Index \\"
local col2 " & (1) & (2) & (3) & (4) & (5) & (6) \\ \midrule "
local rules \cmidrule(lr){2-3} \cmidrule(lr){5-6}
local est_consume mgroups("Stock Remaining" "Consumption" "Daily Protein (g)" "Cons.", pattern(1 0 1 1 0 1) $mgset  end( `rules')) posthead(`col1' `col2')
esttab using "$out_tables/Table_S8.tex", keep($B_main) $estset $mainstats $ecell `est_consume'
	
* Make graph for picture
if `pics'==1 {
coefplot `vbls' eat_day_pfem_win5 ///
	, keep($B_coef ) xline(0) xline(1) $grset $xset name(gconsume) legend(cols(3)) ///
	p1(label("Stock (Kgs.)") msymbol(o) color(edkblue) ciopts(lcolor(edkblue) lwidth(thin))) ///
	p2(label("Stock (Months)") msymbol(d) color(maroon) ciopts(lcolor(maroon) lwidth(thin))) ///
	p3(label("Pulses Consumed") msymbol(t) color(eltgreen) ciopts(lcolor(eltgreen) lwidth(thin))) ///
	p4(label("Total Protein") msymbol(s) color(brown) ciopts(lcolor(brown) lwidth(thin))) ///
	p5(label("Female Protein") msymbol(X) color(gs0) ciopts(lcolor(gs0) lwidth(thin)))
graph export "$out_figures/Parts/fig_coef_consumption_`pic'.png", replace

graph combine gprod gprofit gconsume, cols(1) graphregion(color(white)) ysize(8)
graph export "$out_figures/Figure_`fno'.png", replace
}
}


*****
* Appendix Regressions
*****
		
if `appx'==1 {
local tabs adopt yield prod earn consume
local hets bl_group bl_pulse_prev bl_wealth_high bl_land_high
local num = 0
local rnum = 9

local col1 " & Main & Farmer Group & Grew Pulses & Above Median & Above Median \\"
local col2 " & Result & Member & Pre-Study & Asset Index & Land Owned \\"
local col3 " & (1) & (2) & (3) & (4) & (5)\\ \midrule "
local rules \cmidrule(lr){3-6}
local est_het mgroups("" "Interaction ($ I=1$)", pattern(0 1 0 0 0) $mgset end( `rules')) posthead(`col1' `col2' `col3')

foreach t in `tabs' {
	snapshot restore `++num'
	
	if `dlasso' {
		* Generate round-specific baseline interactions
		// hh_size i.caste i.block i.m_education bl_pulse_prev bl_wealth_win5
		foreach vbl in $controls {
			if ("`vbl'"=="i.block") continue
			
			* Parse out factor variables
			if substr("`vbl'",2,1)=="." {
				local vname = substr("`vbl'",3,.)
				quietly tab `vname', gen(`vname'_)
				local items = "`vname'_*"
			}
			else local items = "`vbl'"
			
			foreach item of varlist `items' {
			forvalues yr = 1/3 {
				gen dl_y`yr'_`item' = `item'*(year==`yr')
			}
			}
		}
		
		* Double Lasso Controls	
		eststo clear
		foreach vbl of varlist `v_`t'' {
			if ("`vbl'"=="eat_day_pfem_win5") continue
			
			local RV $Rvars
			if (`:list posof "`vbl'" in RV') local samp ever
			else local samp srvy
			
			eststo: dsregress `vbl' $B_main i.block if drop_`samp'==0, controls(dl_*) vce(cluster village)
			di "Outcome: `vbl'"
			di "Controls: `e(controls_sel)'"
			if ("`vbl'"=="eat_day_pfem_win5") quietly sum eat_day_pfem_win5 if year==2 & treat_input==0
			else quietly sum `vbl' if year==1 & treat_input==0
			if (substr("`vbl'",1,3)!="pca") estadd scalar ctrl_mean = `r(mean)'
			
		}
		esttab using "$out_tables/Appendix/rob_`t'_dslasso.tex", keep($B_main) $estset $mainstats `est_`t''
		
		drop dl_*
	}
	
	if `robust' {
	*No controls
	eststo clear
	foreach vbl of varlist `v_`t'' {
		local RV $Rvars
		if (`:list posof "`vbl'" in RV') local samp ever
		else local samp srvy
	
		eststo: reg `vbl' $B_main i.block if drop_`samp'==0, vce(cluster village)
		if ("`vbl'"=="eat_day_pfem_win5") quietly sum eat_day_pfem_win5 if year==2 & treat_input==0
		else quietly sum `vbl' if year==1 & treat_input==0
		if (substr("`vbl'",1,3)!="pca") estadd scalar ctrl_mean = `r(mean)'
	}
	esttab using "$out_tables/Appendix/rob_`t'_uncontrol.tex", keep($B_main) $estset $mainstats `est_`t''

	* Unbalanced panel
	eststo clear
	foreach vbl of varlist `v_`t'' {
		eststo: reg `vbl' $B_main $controls, vce(cluster village)
		if ("`vbl'"=="eat_day_pfem_win5") quietly sum eat_day_pfem_win5 if year==2 & treat_input==0
		else quietly sum `vbl' if year==1 & treat_input==0
		if (substr("`vbl'",1,3)!="pca") estadd scalar ctrl_mean = `r(mean)'
	}
	esttab using "$out_tables/Appendix/rob_`t'_unbalanced.tex", keep($B_main) $estset $mainstats `est_`t''

	* Drop Samastipur
	eststo clear
	foreach vbl of varlist `v_`t'' {
		local RV $Rvars
		if (`:list posof "`vbl'" in RV') local samp ever
		else local samp srvy
	
		eststo: reg `vbl' $B_main $controls if drop_`samp'==0 & district!=2, vce(cluster village)
		if ("`vbl'"=="eat_day_pfem_win5") quietly sum eat_day_pfem_win5 if year==2 & treat_input==0 & district!=2
		else quietly sum `vbl' if year==1 & treat_input==0 & district!=2
		if (substr("`vbl'",1,3)!="pca") estadd scalar ctrl_mean = `r(mean)'
	}
	esttab using "$out_tables/Appendix/rob_`t'_noAKRSP.tex", keep($B_main) $estset $mainstats `est_`t''

	* Break down by IDi treatment
	eststo clear
	foreach vbl of varlist `v_`t'' {
		local RV $Rvars
		if (`:list posof "`vbl'" in RV') local samp ever
		else local samp srvy
	
		eststo: reg `vbl' $B_app $controls if drop_`samp'==0, vce(cluster village)
		if ("`vbl'"=="eat_day_pfem_win5") quietly sum eat_day_pfem_win5 if year==2 & treat_input==0
		else quietly sum `vbl' if year==1 & treat_input==0
		if (substr("`vbl'",1,3)!="pca") estadd scalar ctrl_mean = `r(mean)'
	}
	esttab using "$out_tables/Appendix/rob_`t'_idi.tex", keep($B_app) $estset $mainstats `est_`t''
	
	} // if robust
	
	if `rob_pic' {
	* Do robustness as figures instead of tables
	graph drop _all
	
	foreach vbl of varlist `v_`t'' {
		eststo clear
		
		local RV $Rvars
		if (`:list posof "`vbl'" in RV') local samp ever
		else local samp srvy
		
		if ("`vbl'"=="pulse_time_prev") local regs "B_y1_t B_y2_t B_y2"
		else if ("`vbl'"=="eat_day_pfem_win5") local regs "B_y2_t B_y3_t B_y3"
		else local regs "$B_main"
		
		// Run robusness regressions
		eststo `vbl'_m: reg `vbl' $B_main $controls if drop_`samp'==0, vce(cluster village)
		eststo `vbl'_u: reg `vbl' $B_main $controls, vce(cluster village)
		eststo `vbl'_s: reg `vbl' $B_main $controls if drop_`samp'==0 & district!=2, vce(cluster village)
		local extra C_y1_t C_y2_t C_y3_t
		rename (B_*_t B_*_e B_*_s) (hold_*_t B_*_t C_*_t)
		eststo `vbl'_h: reg `vbl' $B_main `extra' $controls if drop_`samp'==0, vce(cluster village)
		rename (B_*_t C_*_t) (C_*_t B_*_t)
		eststo `vbl'_l: reg `vbl' $B_main `extra' $controls if drop_`samp'==0, vce(cluster village)
		rename (hold_*_t C_*_t B_*_t) (B_*_t B_*_e B_*_s)
		
		coefplot `vbl'_m `vbl'_s `vbl'_h `vbl'_l `vbl'_u ///
			, keep($B_coef ) order($B_coef ) xline(0) $grset name(`vbl') title("`: variable label `vbl''") ///
			p1(label("Main Specification") msymbol(X) color(gs0) ciopts(lcolor(gs0) lwidth(thin))) ///
			p2(label("Drop Samastipur") msymbol(o) color(edkblue) ciopts(lcolor(edkblue) lwidth(thin))) ///
			p3(label("Yr. 1 Extension") msymbol(d) color(maroon) ciopts(lcolor(maroon) lwidth(thin))) ///
			p4(label("Yr. 1 Subsidy Only") msymbol(t) color(eltgreen) ciopts(lcolor(eltgreen) lwidth(thin))) ///
			p5(label("Unbalanced Panel") msymbol(s) color(brown) ciopts(lcolor(brown) lwidth(thin)))
			
		graph export "$out_figures/Parts/fig_robust_`vbl'.png", replace
	}
	
	grc1leg `v_`t'', cols(2) graphregion(color(white)) ysize(8)
	graph export "$out_figures/Figure_S`rnum++'.png", replace
	graph drop _all
	} // if rob_pic
			
	*** Heterogeneity interaction tables
	if `hetero' {
	foreach vbl of varlist `v_`t'' {
		local RV $Rvars
		if (`:list posof "`vbl'" in RV') local samp ever
		else local samp srvy
	
		
		eststo clear
		eststo: reg `vbl' $B_main $controls if drop_`samp'==0, vce(cluster village)
		if ("`vbl'"=="eat_day_pfem_win5") quietly sum eat_day_pfem_win5 if year==2 & treat_input==0
		else quietly sum `vbl' if year==1 & treat_input==0
		if (substr("`vbl'",1,3)!="pca") estadd scalar ctrl_mean = `r(mean)'
		
		foreach het in `hets' {
						
			* Generate heterogeneity interactions
			cap drop *_i
			gen byte int_i = `het'
			foreach coef in $B_main {
				gen byte `coef'_i = `coef'*`het'
			}			
			
			* Label variables
			label var int_i "$ I=1$
			label var B_y2_i "Y2*($ I=1$)"
			label var B_y3_i "Y3*($ I=1$)"
			label var B_y1_t_i "Treat Y1*($ I=1$)"
			label var B_y2_t_i "Treat Y2*($ I=1$)"
			label var B_y3_t_i "Treat Y3*($ I=1$)"
			
			* Run regressions
			eststo: reg `vbl' $B_het $controls if drop_`samp'==0, vce(cluster village)
			if ("`vbl'"=="eat_day_pfem_win5") {
				quietly sum `vbl' if year==2 & treat_input==0 & `het'==0
				estadd scalar ctrl_mean = `r(mean)'
				if ("`het'"!="bl_group") {
					quietly sum `vbl' if year==2 & treat_input==0 & `het'==1
					estadd scalar het_mean = `r(mean)'
				}
				else estadd scalar het_mean = .
			}
			else {
				quietly sum `vbl' if year==1 & treat_input==0 & `het'==0
				if (substr("`vbl'",1,3)!="pca") estadd scalar ctrl_mean = `r(mean)'
				if ("`het'"!="bl_group") {
					quietly sum `vbl' if year==1 & treat_input==0 & `het'==1
					if (substr("`vbl'",1,3)!="pca") estadd scalar het_mean = `r(mean)'
				}
				else estadd scalar het_mean = .
			}
			
			* Reset variables
			// rename `het' `vbl'
			// rename interact `het'
		}
		
		esttab using "$out_tables/Appendix/het_`vbl'.tex", keep($B_het) order($B_het) $estset $mainstats `est_het'
		// $hetstats
	} // foreach vbl
	} // if hetero
	
}
}
snapshot erase _all



